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INTRODUCTION 


Currently,  the  Naval  Air  Warfare  Center  Weapons  Division  (NAWCWPNS)  is 
exploring  using  the  solution  of  the  time-dependent  Schrodinger  for  some  applications  in 
signal  processing.  The  approach  being  taken  to  solve  this  problem  is  a  difference  equation. 
However,  the  solution  based  upon  a  difference  equation  is  unstable,  and  the  solution  must 
be  renormalized  every  time  step.  This  report  reviews  two  methods  to  obtain  solutions  that 
are  stable  and  unitary  (preserve  the  norm).  Both  methods  are  based  upon  the  split-operator 
approach.  One  meAod,  called  the  k-space  method,  will  use  the  fast  Fourier  transform 
(FFT),  and  one  method,  called  the  R-space  method,  will  not.  The  A:-space  method  is  more 
accurate  than  the  R-space  method.  However,  it  is  a  “global”  method  (because  it  uses  the 
FFT),  while  the  i?-space  method  is  “local.”  This  difference  can  be  exploited  to  efficiently 
use  tile  available  computing  architecture.  Included  in  this  report  are  two  software  packages 
(one  in  C  and  one  in  C++)  that  implement  these  methods.  The  outline  of  this  report  is  as 
follows:  section  I  describes  the  methods  for  calculating  the  time  evolution  of  the  wave 
function  governed  by  a  time  dependent  Schrodinger  equation,  the  Potential  Section 
describes  the  potentials  that  are  used  in  the  Schr(^nger  equation,  and  the  Software 
Packages  Section  contains  the  two  software  listings. 


THE  TIME-DEPENDENT  SCHRODINGER  EQUATION 

The  equation  that  governs  the  time  evolution  of  our  wave  function  is  the  time-dependent 
Schrodinger  equation 


2  2 

^  ^  \i/ix,t)  +  V(x,t)ii/(xj)  (1) 

ot  2m 


where  V(x,t)  is  the  potential.  Given  the  initial  wave  function  we  can  get  the  evolution  of 
the  wave  function  in  time  by  integrating  Equation  1  over  a  small  time  step  where  the 
Hamiltonian  H  is  approximately  constant.  In  such  case  the  solution  is  straightforward 

y/(x,t+At)  =  (2) 


This  is  repeated  to  get  the  time  evolution  for  all  time.  Note  that  Equation  2  is  unitary  if  the 
Hamiltonian  is  real,  i.e.  the  probability  is  conserved,  and  we  do  not  need  to  renormalize  the 
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wave  function  after  each  time  step.  The  right  hand  side  of  Equation  2  is  evaluated  as 
follows. 


This  can  be  seen  by  expanding  the  exponential  on  both  sides  of  Equation  3  and  compare 
term  by  term.  Equation  3  is  the  essence  of  the  split-operator  approach. 


THE  k-SPACE  METHOD 

The  first  and  last  exponential  in  Equation  3  are  diagonal  in  real  space  and  can  be 
evaluated  in  a  straightforward  manner.  The  middle  exponential,  which  contains  the 
Laplacian,  can  be  evaluated  accurately  and  most  efficiently  in  k  space.  Therefore  the  steps 
for  evaluating  Equation  3  are  as  follows: 

a.  Multiply  the  first  exponential  with  the  wave  function, 

b.  Fourier  transform  the  resulting  product  to  ^-space  using  FFT, 

c.  Multiply  with  the  middle  exponential, 

d.  Fourier  transform  back  to  real  space, 

e.  Multiply  with  the  last  exponential. 

Note  that  the  use  of  the  Fourier  transform  on  a  finite  grid  means  that  we  are  using  periodic 
boundary  condition. 


THE  R-SPACE  METHOD 

The  disadvantage  of  the  k-space  method  is  that  the  FFT  is  a  global  operation  and 
therefore  cannot  used  to  exploit  parallel  computing  architectures  (particularly  a 
distributed  environment).  The  R-space  method  (Reference  1)  is  intended  to  remedy  this 
deficiency.  To  evaluate  the  middle  exponential  operator  we  use  the  same  split-operator 

trick.  Let  us  denote  the  Laplacian  operator  as  (L  /Ax^ ),  which  can  be  written  as  a  sum  of  2 
parts:  A  and  5.  We  then  have 


^ihLAt  /(2mAx^)  _  l(4mAx^)^iWAt /(2TnAx^)gitiAAt 


To  understand  how  and  why  we  do  this  it  is  easier  to  look  at  a  simple  example  of  sn  8- 
point  grid.  In  this  case  the  Laplacian  operator,  using  a  three-point  formula  and  periodic 
boundary  condition,  has  the  form 
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The  Laplacian  L  can  be  written  as  the  sum  of  two  matrices  A  and  B  as  follows 
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Notice  that  A  and  5  are  2  x  2  block  diagonal;  therefore,  the  exponentiation  of  these 
operators  can  be  evaluated  easily.  The  2x2  matrix  that  we  need  to  diagonalize  is 


M  = 


(7) 


This  matrix  has  eigenvalues  of  0  and  -2.  The  matrix  that  diagonalize  M  is 


S  =  S~^ 


n 


(8) 


That  is 
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SMS~^  = 


^0 

.0 


(9) 


The  block  diagonal  nature  of  the  matrices  makes  Equation  4  very  easy  to  evaluate.  The 
problem  of  diagonalizing  a  NxN  matrix  is  reduced  to  the  problem  of  diagonalizing  (N/2) 
identical  2x2  matrices  (which  we  already  done  in  Equations  7  through  9).  This  algorithm 
is  highly  parallel  as  the  2  x  2  blocks  can  be  done  in  parallel.  For  each  2x2  block  we  have 


(10) 


A  final  note  about  the  R-space  method  is  that  it  can  be  generalized  to  higher  dimension  as 
the  Laplacian  operator  for  each  dimension  commutes  among  themselves  (i.e.  =  L^L^^). 


THE  POTENTIAL 


THE  BASIC  POTENTIAL 

In  the  Time-Dependent  Schrodinger  Equation  Section,  we  discussed  the  solution  to  the 
Schrodinger  equation  for  an  arbitrary  potential.  In  this  section  we  briefly  describe  the  types 
of  potential  available  in  the  software  packages.  This  section  is  intended  to  provide  some 
scientific  background  to  the  problem.  The  programs  allow  the  specification  of  1  or  2 
potentials.  The  parameters  to  be  specified  for  each  potential  are: 

1.  Its  center 

2.  Its  width 

3.  Its  velocity 

4.  Its  depth  (height) 

5.  Its  shape:  Square  or  l/(cosh^) 

Further  description  of  these  parameters  can  be  found  in  the  Software  Packages  Section. 
Here  we  want  to  give  some  details  regarding  item  5.  The  use  of  a  1/cosh^  potential  will 
help  us  verify  the  program  as  well  as  gain  insight  into  the  physical  process  (for  actual 
application)  since  the  bound  states  of  this  potential  are  known.  This  potential  has  the  form 
V(x)  =  -V^cosh^(xlw).  We  will  just  state  the  results  for  the  bound  states  here;  readers 
interested  in  the  full  solution  should  read  Reference  2.  The  number  of  bound  states  for  this 
potential  is  n+ 1  where  n  is  the  largest  integer  smaller  than  s  with 


2 


1  + 


fW 


(11) 
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The  bound  states  are 


Ynix)  =  (l - ^  )  F(-n,2s  -n  +  l,s-n  + 1;(1  -  ^)/2) 


(12) 


where  ^  =  tanh(x/w), «  =  0, 1, 2, and  F  is  the  hypergeometric  function 


Ftofcc;.)  =  1  ^ ^  ^  ^ 

c  2c(c  +  l)  2-3-c(c  +  lXc  +  2) 


The  energy  of  the  nth  bound  state  is 


E(n)  =  -- - 

Imw^ 


(14) 


From  Equation  1 1  we  can  choose  the  potential  parameters,  w  and  V^,  and  the  mass  m  so 
that  the  well  has  only  one  bound  state.  Since  the  mass  m  is  arbitrary  we  define  it  in  terms 

of  a  new  variable  j3, 


"  '^max  _  y 


2m 


o- 


(15) 


The  variable  P  is  the  ratio  of  the  maximum  potential  energy  to  the  maximum  kinetic  energy. 
For  a  grid  with  N  points  and  length  L,  =  (tiNIL). 

The  program  has  an  option  in  which  one  can  specify  the  initial  wave  function  to  be  a 
constant  or  die  lowest  bound  state  (as  calculated  from  Equation  12).  When  the  bound  state 
is  chosen  with  appropriate  setting  of  the  potential  parameters  (one  potential  with  a  1/cosh^ 
shape  and  zero  velocity),  the  wave  function  should  remain  in  this  state  for  all  time.  This 
can  be  verified  by  computing  the  overlap  of  the  wave  function  with  this  bound  state  as  a 
function  of  time. 

In  addition  to  the  real  potential,  the  program  allows  the  addition  of  an  absorptive 
(imaginary)  potential  at  both  ends  of  the  grid.  The  purpose  of  this  potential  is  to  cause  die 
wave  function  to  go  to  zero  at  the  boundary  of  the  grid  and  prevent  aliasing  or  wrap-around 
effects  of  the  FFT.  Unfortunately,  this  potential  also  scatters  some  of  the  wave  function 
back.  An  ideal  situation  is  to  use  some  transmitting  boundary  condition  without  reflection, 
but  this  is  not  known  for  the  Schrodinger  equation  at  present. 
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HILDE  FILTER 

Besides  the  basic  potential  one  can  do  the  following  operations  on  the  potential  before 
using  it  in  the  Schrodmger  equation: 

1.  Add  random  noise. 

2.  Pass  it  through  a  Hilde  filter  (developed  by  Jeff  Hilde,  NAWCWPNS),  to  be 
rectified,  amplified,  and  convolved  with  a  Gaussian. 

Figure  1  shows  the  flow  diagram  for  the  whole  program,  and  Figure  2  shows  the  steps  to 
generate  the  potential  for  use  in  the  Schrodinger  equation. 


FIGURE  1.  Flow  Chart  of  the  Program. 
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FIGURE  2.  Flow  Chart  for  the  Generation  of  the  Potential  for  the  Schrodinger  Equation. 


SOFTWARE  PACKAGES 


This  section  contains  two  source  codes:  one  in  C  and  one  in  C++.  The  code  in  C  was 
developed  with  Borland  C++  (4.5)  compiler  running  on  a  486  PC.  If  you  are  not 
developing  your  own  code  then  Ais  code  is  adequate.  All  the  changeable  parameters  (with 
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comments/descriptions  in  bold)  are  at  the  beginning  of  the  program.  This  program  will  give 
a  graphical  output  showing  the  time  evolution  of  the  potential  and  the  wave  ftmction. 

If  you  are  developing  your  own  code  to  solve  the  Schrodinger  equation  for  a  potential 
of  your  choice  then  all  you  want  is  the  time  propagation  function.  This  function,  in  its 
most  useful  form,  is  found  in  the  second  program  written  in  C++.  To  demonstrate  how  it 
works,  it  was  wrapped  within  a  program  similar  to  the  one  written  in  C  (but  without  the 
graphical  outputs).  This  program  was  developed  with  a  C++  compiler  running  on  an  SGI. 
The  time  propagation  function  has  3  parts:  the  FFT  function  fourl,  the  function 
Smultiply,  and  the  function  time_stepping.  You  can  substitute  your  own  FFT  routine 
if  you  desire.  These  three  functions  wiU  solve  the  Schrodinger  equation  using  either  the  k- 
space  orR-space  method. 

Both  programs  can  be  obtained  from  the  authors  if  you  do  not  feel  like  retyping  the 
program. 


C  PROGRAM 

#include  <graphics.h> 

#include  <stdlib.h> 

#include  <stdio.h> 

#include  <conio.h> 
tinclude  <math.h> 

# include  <dos.h> 

#include  <fstream.h> 

/  /  for  random  number  generator 

#define  Ml  259200 
#define  lAl  7141 
#define  ICl  54773 
#define  RMl  (1.0 /Ml) 

#define  M2  134456 
#define  IA2  8121 
#define  IC2  28411 
#define  RM2  (1.0/M2) 

#definG  M3  243000 
idefine  IA3  4561 
#define  IC3  51349 

//  global  variables 

#define  nsize  128 
#define  Hfactor  0.8 

float  v[nsi2e],  v0[nsi2e],  hil[nsi2e],  hrec[nsize],  vnoise [nsize] ; 
float  psilr,  psili,  psi2r,  psi2i; 

float  k2rQ[nsize]  ,  potr[nsize],  poti[nsize],  dat [2*nsi2e] ; 
float  ground_state [nsize]  ; 

float  pi  =  3.141592653589  ; 


//  Number  of  grid  points  (pixels) 

//  This  is  to  tune  the  Hilde  filter 
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float  gain  =  80 . 
float  dt  =  .05  ; 
int  nt  =20  ; 
long  seed  =  3411 


//  For  hllde  filter 
//  Size  of  time  step 

//  Number  of  time  step  between  frame 
//  Seed  for  random  number  generator 


/*  Whenever 

the  mass 

is 

changed  (make  smaller) 

remember  to 

check  that 

the  time 

step  dt 

is  adequate  !  !  !  M  1 

!  i  */ 

float  beta  = 

o 

o 

// 

Mass  parameter 

float  length 

=  12.8  ; 

// 

Length  of  grid 

//  Options 

int  nopt  =  0 

/ 

// 

1  = 

1  potential,  0=2 

potentials 

int  iopt  =  0 

/ 

// 

1  = 

square,  0  =  1/cosh 

^2  potential 

int  potopt  = 

0  ; 

// 

=  1 

use  hilde  filter. 

=  0  don '  t 

int  wopt  =  0 

/ 

// 

1  = 

const.  psi(t=0),  0 

=  ground  state 

int  output  = 

0  ; 

// 

Output  overlap  data  to 

a  file ;  1  =  yes 

int  kspace  = 

1  ; 

// 

Using  k-space  (=1)  or 

R-space  (=0) 

//  Data  about  potentials 


float 

vel  =  0  ; 

// 

Vel.  of 

potential 

(pixel /frame) 

float 

xcenter  =  55  ; 

// 

Center 

of 

potential 

initially 

float 

vdepth  =  0.5  ; 

// 

Depth 

of 

potential 

float 

width  =  4  ; 

// 

Width 

of 

potential 

in  pixels 

float 

vel2  =  0  ; 

// 

Velocity 

of  2nd  potential 

float 

xcenter2  =70  ; 

// 

Center 

of 

2nd  potential  initially 

float 

vdepth2  =1.0  ; 

// 

Depth 

of 

2nd  potential 

float 

width2  =  4  ; 

// 

Width 

of 

2nd  potential  in  pixels 

float 

wa  =  4  ; 

// 

Width 

of 

absorbing 

potential 

float 

Va  =  10  ; 

// 

Depth 

of 

potential. 

0  =  no  pot. 

float 

noise_factor  =  1  ; 

// 

0  =  no 

>  noise  ,  1  = 

yes 

float 

snr  =  5  ; 

// 

Signal 

to 

noise 

//x-y  coordinate  and 

int  x_target  =30  ; 
int  y_target  =10  ; 
float  scale^target  =  5  ; 
int  x_noise  =30  ; 
int  y__noise  =  100  ; 
float  scale_noise  =  5  ; 
int  x_hilde  =  30; 
int  y__hilde  =  200; 
float  scale_hilde  =  100; 
int  x_jpot  =30  ; 
int  y__pot  =  250; 
float  scale_pot  =  5  ; 
int  x_prob  =  310; 
int  y_prob  =  200; 
float  scale_prob  =  100  ; 


scale  of  the  plots 
//  Position  of 

//  Position  of 

//  Pos . 

//  Pos . 

//  Position 


potential 


(potential  +noise) 


of  output  from  Hilde  filter 


of  the  Schrodinger  potential 


of  probability 
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//Gaussian  for  convolution 

float  gauss[9]  =  {54,129,242,352,399,352,242,129,54}  ; 

//  Cone  Filter 

float  cone  [9]  =  {l.,2.,3.,4.,5.,4.,3-,2,,l.}  ; 

//  Random  number  generator  (Adapted  from  "Numerical 
Recipes... by  Press  et.  al.") 

float  rani (long  idum) 

{ 

static  long  ixl,  ix2,  ix3  ; 
static  float  r[98]  ; 
float  temp  ; 
static  int  iff=0  ; 
int  j  ; 

if  (idum  <0  M  iff  ==  0) 

{ 

iff  =  1  ; 

ixl  =  (ICl  -  idum)  %  Ml  ; 

ixl  =  (IAl*ixl  +  ICl)  %  Ml  ; 

ix2  =  ixl  %  M2  ; 

ixl  =  (IAl*ixl  +  ICl)  %  Ml  ; 

ix3  =  ixl  %  M3  ; 

for  (j=:l;  j<=97;  j++) 

{ 

ixl  =  (IAl*ixl  +  ICl)  %  Ml  ; 
ix2  =  (IA2*ix2  +  IC2)  %  M2  ; 
r[j]  =  (ixl  +  ix2*RM2)*RMl  ; 

} 

idum  =  1  ; 

} 

ixl  =  (IAl*ixl  +  ICl)  %  Ml  ; 

ix2  =  (IA2*ix2  +  IC2)  %  M2  ; 

ix3  =  (IA3*ix3  +  IC3)  %  M3  ; 

j  =  1  +  (  (97*ix3)/M3)  ; 

temp  =  r[j]  ; 

r[j]  =  (ixl+ix2*RM2)  *RM1  ; 
return  temp  ; 


//  Graphics  function 

void  Draw(int  x,  int  y,  float  d[nsize]  ,  float  scale) 

{ 

int  i,j; 
j  =  x; 

moveto  (x,y~d[0]  *5*scale)  ; 
for  (i=l;  i<nsize;  i++) 

{ 

j  =  j+  (512  /  nsize/2); 
lineto(  j  ,y’"d[i]  *5*scale)  ; 
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} 

} 

//  Convolve  input  f  (x)  with  a  filter  g(x)  to  get  h(x) 

void  Conv(float  f[nsi2e],  float  g[],  float  h[nsize] ,  int  gsize) 

{ 

int  i,  j  ; 
float  sum  ; 

for  (i=0;  i<nsize;  i++) 

{ 

sum  =  0 .  ; 

for  ( j=- (gsize-1) ;  j  <  gsize;  j++) 

{ 

sum  =  sum  +  g[j+gsize-l]  *  f [ {nsize+i+ j ) %nsize]  ; 

} 

h  [  i  ]  =  sum  ; 

} 

} 

/  /  Hilde  filter 

void  hilde (float  xl[nsize],  float  xO[nsize],  float  out[nsi2e],  int  flag) 

{ 

static  float  diff[nsize],  al[nsize],  a2[nsize],  a3 [nsize]  ; 
int  i  ; 

if  (flag  ==  1) 

{ 

for  (i=0;  i<nsize;  i++) 

{ 

al[i]  =  0  ; 
a2[i]  =  0  ; 
a3[i]  =  0  ; 

} 

} 

for  (i=0;  i<  nsize;  i++) 

{ 

diff[i]  =  xl[i]  -  xO[i]  ; 

} 

//  convolution  with  a  cone  filter 

Conv(dif f , cone, out , 5)  ; 

for  (i=0;  i<nsize;  i++) 

{ 

al[i]  =  Hfactor*al [i]  +  ( 1-Hf actor) *out [ i]  ; 
a2[i]  =  Hfactor*a2 [i]  +  ( l-Hfactor) *al [ i ]  ; 
a3[i]  =  Hfactor*a3 [i]  +  ( 1-Hfactor) *a2 [ i ]  ; 
out [ i ]  =  a3 [ i ]  ; 

xO [i]  =  xl[i]  ; 

} 
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} 


//  rectify  the  output  from  the  Hilde  filter 

void  rect( float  in[nsize],  float  out[nsize]) 

{ 

int  i  ; 

float  tem[nsize]  ; 

//  Rectification 

for  (i=0  ;  i<nsi2e;  i++) 

{ 

if  (in[i]  <  0) 

tem[i]  =  gain*in[i]  ; 

else 

temfi]  =  0  ; 

} 


} 


//  convolve  with  a 

Conv(tein,  gauss,  out ,  5) 
Conv ( out , gaus  s , t  em, 5 ) 
Conv  ( tern,  gauss ,  out ,  5 ) 


gaussian  3  times 


//  fft  routine  (Adapted  from  "Numerical  Recipes*.. by  Press  et . 
al .  ") 

void  fourKfloat  data[],  int  nn,  int  isign) 

{ 

int  n,  mmax,  m,  j,  istep,  i  ; 

double  wtemp,  wr,  wpr,  vzpi,  wi,  theta  ; 

float  tempr,  tempi,  tern  ; 


n  = 

j  = 

for 


{ 


nn  «  1  ;  //n  =  nn*2 
1  ; 

(i=l;  i<n;  i+=2) 


if  (j  >  i) 

{ 

tern  =  data[i-l]  ; 
data[i-l]  =  data[j“l]  ; 
data[j*-l]  =  tern  ; 
tern  =  data[i]  ; 
data[i]  =  data[j]  ; 
data[j]  =  tern  ; 

} 

m  =  n  »  1  ;  //m  =  n/2 
while  (m>=2  ScSc  j  >  m) 

{ 

j  -=  m  ;  //j  =  j  -  m 

m  »=  1  ;  //m  =  m/2 

} 

j  +=  m  ;  //j  =  j+m 
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} 

mmax  =  2  ; 
while  (n  >  mmax) 

{ 

i s  t  ep  =  2  *mmax  ; 

theta  =  6.28318530717959/ (isign*mmax)  ; 

wtemp  =  sin ( 0 . 5*theta)  ; 

wpr  =  -2 . 0*wtemp*wtemp  ; 

wpi  =  sin (theta)  ; 

wr  =  1.0  ; 

wi  =  0.0  ; 

for  (m=l;  m  <  mmax;  m  +=  2) 

{ 

for  (i=m;  i<=n;  i  +=  istep) 

{ 

j  =  i  +  mmax  ; 

tempr  =  wr*data[j”l]  “  wi*data[j] 
tempi  =  wr*data[j]  +  wi*data[j“l] 
data[j-l]  =  data[i-'l]  -  tempr  ; 
data[j]  =  data[i]  -  tempi  ; 
data[i-l]  +=  tempr  ; 
data[i]  +=  tempi  ; 

} 

wr  =  (wtemp=wr) *wpr  -  wi*wpi  +  wr  ; 
wi  =  wi*wpr  +  wtemp*wpi  +  wi  ; 

} 

mmax  =  istep  ; 

} 

//if  isign  =  -1  multiply  by  1/nn 
if  (isign  ==  -1) 

{ 

for  (i=0;  i  <  2*nn;  i++) 

{ 

data[i]  =  data[i]/nn  ; 

} 

} 

} 

//  function  to  multiply  with  the  2x2  S-matrix 

void  SmultiplyO 

{ 

double  templr,  tempi i,  temp2r,  temp2i  ; 
double  sll  =  l./sqrt(2.); 
double  sl2  =  l./sqrt(2.); 

double  s21  =  l./sqrt(2.); 

double  s22  =  -1 . /sqrt (2 . ) ; 

templr  =  sll*psilr  +  sl2*psi2r  ; 
templi  =  sll*psili  +  sl2*psi2i  ; 
temp2r  =  s21*psilr  +  s22*psi2r  ; 
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teinp2i  =  s21*psili  +  s22*psi2i  ; 
psilr  =  (float) (templr)  ; 

psi2r  =  (float) (temp2r)  ; 

psili  =  (float) (tempi!)  ; 

psi2i  =  (float) (temp2i)  ; 


//  Time  propagation  function 

void  time_stepping ( float  psir[nsize],  float  psii[nsize],  int  nn,  float 
mass) 

{ 

int  i,  j  ; 

float  alpha,  temr,  temi,  dx,  tem2r,  tem2i  ; 
dx  =  length/nsize  ; 

//  applied  exp[~iVt/2]  to  the  wave  function 

for  (i=0;  i<nn  ;  i++) 

{ 

alpha  =  potr[i]*dt/2  ; 

temr  =  cos (alpha) *psir [ i]  +  sin (alpha) *psii [ i ]  ; 

temi  =  cos (alpha) *psii [ i]  -  sin (alpha) *psir [ i]  ; 
psir[i]  =  temr*exp(“poti [i] *dt/2)  ; 
psii[i]  =  temi*exp(-poti [i] *dt/2)  ; 

} 

if  (kspace) 

{ 

//  first  fft  to  k  space 

for  (i=0;  i<nn  ;  i++) 

{ 

j  =  2*(i)  ; 

dat[j]  =  psir[i]  ; 
dat[j+l]  =  psii[i]  ; 

} 

fourl (dat , nn, 1)  ; 
for  (i=0;  i  <nn  ;  i++) 

{ 

j  =  2*(i)  ; 
psir[i]  =  dat[j]  ; 
psii[i]  =  dat[j+l]  ; 

} 

//  applied  exp [-i (k**2) t/2m]  to  the  wave  function 

for  (i=0;  i<nn  ;  i++) 

{ 

alpha  =  k2m[i]*dt/2  ; 

temr  =  cos (alpha) *psir[i]  +  sin (alpha) *psii [ i]  ; 
temi  =  cos (alpha) *psii [ i]  -  sin (alpha) *psir[i]  ; 
psir[i]  =  temr  ; 
psiiti]  =  temi  ; 
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} 

//  fft  back  to  real  space 

for  (i=0;  i<nn  ;  i++) 

{ 

j  =  2*(i)  ; 
dat[j]  =  psir[i]  ; 
dat[j+l]  =  psii[i]  ; 

} 

fourl (dat, nn, -1)  ; 
for  (i=:0;  i  <nn  ;  i++) 

{ 

j  =  2*(i)  ; 
psir[i]  =  dat[j]  ; 
psii[i]  =  dat[j+l]  ; 

} 

} 

else 

{ 

//  everything  is  done  in  real  space 
//  operate  exp [iAdt / (4*m*dx*dx) ] 

for  {i=0;  i<nsize  ;  i  +=  2) 

{ 

psilr  =  psir[i]  ; 
psili  =  psii[i]  ; 
psi2r  =  psir[i+l]  ; 
psi2i  =  psii[i+l]  ; 

SmultiplyO  ; 

alpha  =  dt/4/inass/dx/dx*  (-2)  ; 

tem2r  =  cos (alpha) *psi2r  -  sin (alpha) *psi2i 

tein2i  =  cos  (alpha)  *psi2i  +  sin  (alpha)  *psi2r 

psi2r  =  tem2r  ; 

psi2i  =  tem2i  ; 

SmultiplyO  ; 
psir[i]  =  psilr  ; 
psii  [i]  =  psili  ; 
psir[i+l]  =  psi2r  ; 
psii[i+l]  =  psi2i  ; 

} 

//  operate  exp  [iBdt/ (2*in*dx*dx)  ] 

psilr  =  psir[0]  ; 
psili  =  psii [0]  ; 
psi2r  =  psir [nsize-l]  ; 
psi2i  =  psii [nsize-1]  ; 

SmultiplyO  ; 

alpha  =  dt/2/mass/dx/dx* ( -2 )  ; 

tem2r  =  cos (alpha) *psi2r  -  sin (alpha) *psi2i  ; 

tem2i  =  cos  (alpha)  *psi2i  +  sin  (alpha) '*^psi2r  ; 

psi2r  =  tem2r  ; 

psi2i  =  tem2i  ; 
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SmultiplyO  ; 
psir[0]  =  psilr  ; 
psii[0]  =  psili  ; 
psir [nsize-1]  =  psi2r  ; 
psii [nsize-1]  =  psi2i  ; 

for  (i=l;  i<nsize-l  ;  i  +=  2) 

{ 

psilr  =  psir[i]  ; 
psili  =  psii[i]  ; 
psi2r  =  psir[i+l]  ; 
psi2i  =  psii[i+l]  ; 

SmultiplyO  ; 

alpha  =  dt/2/mass/dx/dx*  (’-2)  ; 

tem2r  =  cos  (alpha)  *psi2r  sin  (alpha)  *psi2i  ; 

tem2i  =  cos (alpha) *psi2i  +  sin (alpha) *psi2r  ; 

psi2r  =  tem2r  ; 

psi2i  =  tem2i  ; 

SmultiplyO  ; 
psir[i]  =  psilr  ; 
psii [i]  =  psili  ; 
psir[i+l]  =  psi2r  ; 
psii[i+l]  =  psi2i  ; 

} 

//  operate  exp  [iAdt  /  {4*in*dx*dx)  ] 

for  (i=0;  i<nsize  ;  i  +=  2) 

{ 

psilr  =  psir[i]  ; 
psili  =  psii[i]  ; 
psi2r  =  psir[i+l]  ; 
psi2i  =  psii[i+l]  ; 

SmultiplyO  ; 

alpha  =  dt/4/mass/dx/dx* (“2)  ; 

tem2r  =  cos (alpha) *psi2r  -  sin (alpha) *psi2i  ; 

tem2i  =  cos (alpha) *psi2i  +  sin (alpha) *psi2r  ; 

psi2r  =  tem2r  ; 

psi2i  =  tem2i  ; 

SmultiplyO  ; 
psir[i]  =  psilr  ; 
psii[i]  =  psili  ; 
psir[i+l]  =  psi2r  ; 
psii[i+l]  =  psi2i  ; 

} 

} 

//  applied  expt-iVt/2]  to  the  wave  function  again 

for  (i=0;  i<nn  ;  i++) 

{ 

alpha  =  potr[i]*dt/2  ; 

temr  =  cos (alpha) *psir [ i]  +  sin (alpha) *psii [i]  ; 
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temi  =  cos  (alpha)  *psii  [  i  ]  -  sin  (alpha)  """psir  [  i  ]  ; 
psir[i]  =  teinr*exp(-poti  [i]  *dt/2)  ; 
psii[i]  =  temi*exp (“poti [i] *dt/2)  ; 

} 

} 

//  function  to  generate  the  potential 

void  potent ial_generat ion (long  int  t) 

{ 

double  xupper,  xlower,  x,  xupper2,  xlower2,  y,  x2,  y2  ; 
int  i  ; 

if  (nopt)  //  option  for  1  or  2  wells 

{ 

if  (iopt)  //  option  for  square  or  l/cosh^2  well 

{ 

//  set  up  square  potential 

xlower  =  vel*t  -  width/2  +  xcenter  ; 
xupper  =  vel*t  +  width/2  +  xcenter  ; 
xlower  =  fmod (( xlower +nsize)  ,  nsize)  ; 
xupper  =  finod(  (xupper+nsize)  ,  nsize)  ; 
for  (i=0  ;  i  <  nsize  ;  i++) 

{ 

v[i]  =  0  ; 

if  (xlower  <  xupper) 

{ 

if  (i  >=  xlower  StSc  i  <=  xupper)  v[i]  =  -vdepth; 

} 

else 

{ 

if  (i  <=  xupper  | |  i  >=  xlower)  v[i]  =  -vdepth; 

} 

} 

} 

else 

{ 

//  set  up  l/cosh**2  potential 

for  (i=0;  i<nsize;  i++) 

{ 

X  =  i  -  xcenter  -  vel*t  ; 
y  =  (cosh(x/width) )  ; 
v[i]  =  -vdepth/y/y  ; 

} 

} 

) 

else 

{ 

if  (iopt)  //  option  for  square  or  l/cosh^2  well 

{ 

//  square  well 

xlower  =  vel*t  -  width/2  +  xcenter  ; 
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xupper  =  vel*t  +  width/ 2  +  xcenter  ; 
xlower  =  fmod ( (xlower+nsize)  ,  nsize)  ; 
xupper  =  fmod( (xupper+nsize)  ,  nsize)  ; 

xlower2  =  vel2*t  -  width2/2  +  xcenter2  ; 
xupper2  =  vel2*t  +  width2/2  +  xcenter2  ; 
xlower2  =  fmod ( (xlower2+nsize)  ,  nsize)  ; 
xupper2  =  fmod ( {xupper2+nsize)  ,  nsize)  ; 

for  (i  =  0  ;  i  <  nsize  ;  i4-+) 

{ 

v[i]  -  0  ; 

//  check  if  point  falls  into  1st  well 

if  (xlower  <  xupper) 

{ 

if  (i  >=  xlower  ScSc  i  <=  xupper)  v[i]  =  -vdepth; 

} 

else 

{ 

if  (i  <=  xupper  |  |  i  >=  xlower)  v[i]  =  -vdepth; 

} 

//  check  to  see  if  point  falls  into  2nd  well 

if  (xlower2  <  xupper2) 

{ 

if  (i  >=  xlower2  ScSc  i  <=  xupper2 ) .  v  [  i  ]  = 

-vdepth2 ; 

} 

else 

{ 

if  (i  <=  xupper 2  M  i  >=  xlower 2)  v[i]  = 

-vdepth2 ; 

} 

} 

} 

else 

{ 

//  set  up  l/cosh**2  potential 

for  (1=0;  i<nsize;  i+-f-) 

{ 

X  =  i  -  xcenter  -  vel*t  ; 

y  =  (cosh(x/width) )  ; 

x2  =  i  -  xcenter 2  -  vel2*t  ; 

y2  =  (cosh (x2/width2 )  )  ; 

v[i]  =  -  vdepth/y/y  -  vdepth2/y2/y2  ; 


//  add  noise  to  potential 
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void  add_noise  (float  snr,  long  idum) 

{ 

float  X,  y  ; 
int  i  ; 

X  =  pow(10 . , snr/10 . )  ; 
y  =  vdepth/x  ; 

for  (i=rO;  i<nsize;  i++) 

( 

vnoise[i]  =  v[i]  +  (rani (idum) -0 . 5) *2*y*noise_factor  ; 

} 

} 

//  Main  program 

void  main  (void) 

{ 

int  i,  flag,  it  ; 
long  int  t; 

float  sum,  dk,  mass,  factl,  fact2,  fr,  fi,  res  ; 
float  temp[nsize],  psir[nsize],  psii[nsize],  prob[nsize]  ; 
double  z,  ss,  yy  ; 
long  idum  ; 

idum  =  seed  ; 

of stream  f out ( "overlap . txt “ )  ; 
int  page=l ; 
int  idelay=0; 

/*  request  autodetection  */ 

int  gdriver  =  EGA,  gmode  =  EGAHI,  errorcode; 

/*  initialize  graphics  mode  */ 

initgraph(Scgdriver,  Scgmode,  "c :  \\bc45\\bgi " )  ; 


/*  read  result  of  initialization  */ 

errorcode  =  graphresult ( ) ; 

if  (errorcode  1=  grOk)  /*  an  error  occurred  */ 

{ 

print f ( "Graphics  error:  %s\n",  grapherrormsg (errorcode) ) ; 
print f ( "Press  any  key  to  halt : " ) ; 
getch ( ) ; 

exit(l);  /*  return  with  error  code  */ 

} 

//  normalize  gaussian 

sum  =  0  ; 
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for  (i=0;  i<9;  i++) 

sum  =  sum  +  gauss [i]  ; 
for  (i=0;  i<9;  i++) 

gauss ti]  =  gauss [i] /sum  ; 

//  normalize  cone  filter 

sum  =  0  ; 

for  (i=0;  i<9 ;  i++) 

sum  =  sum  +  cone [ i ]  ; 

for  (i=0;  i<9;  i++) 

conG[i]  =  cone [i] /sum  ; 

//  calculate  the  mass  from  beta  (a  ratio  of  PE  to  KE) 

mass  =  beta* {pi*nsize/length) * (pi*nsize/length) /2/vdepth  ; 

//  set  up  the  kinetic  energy  vector 

for  {i=l;  i<nsize+l  ;  i++) 

{ 

dk  =  2*pi/length  ; 
k2m[i-l]  =  dk*(i-l)  ; 
if  (i  >  (nsize/2+1)) 

{ 

k2m[i-l]  =  dk* (i-nsize-l)  ; 

} 

k2ra[i-l]  =  k2m[i-l] *k2m[i~l] /mass  ; 

} 

if  (wopt) 

{ 

//  initialize  wave  function  to  a  constant 

for  (i=0  ;  i<nsize  ;  i++) 

{ 

psir[i]  =  1/sqrt ( length)  ; 
psii[i]  =  0  ; 

} 

} 

else 

{ 

//  ground  state  of  l/cosh**2  potential 

for  (i=0;  i<nsize  ;  i++) 

{ 

z  =  tanh (( i-xcenter) /width)  ; 
res  =  nsize/length  ; 

ss  =  0 . 5* (-1+sqrt (l+8*mass*vdepth*width*width/res/res) ) 

yy  =  1  -  z*z  ; 

psir[i]  =  pow(yy, ss/2 . )  ; 

psii[i]  =  0  ; 

} 

sum  =  0 ; 

for  (i=0;  i<nsize  ;  i++) 

{ 
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sum  =  sum  +  psir [ i] *psir [i]  ; 

} 

sum  =  sqrt(sum)  ; 
for  (i=0;  i<nsize;  i++) 

{ 

psir[i]  =  psir [i] /sum  ; 
ground_state [i]  =  psir[i]  ; 

} 


//  set  up  the  absorbing  potential  at  the  end  of  the  grid 

for  (i=0;  i<nsize  ;  i++) 

{ 

factl  =  i*i/wa/wa  ; 

fact2  =  (i-nsize+1) * ( i-nsize+1) /wa/wa  ; 

if  (factl  >  30)  factl  =  30  ; 

if  (fact2  >  30)  fact2  =  30  ; 

poti[i]  =  Va  *  (exp{-factl)  +  exp(-fact2))  ; 

} 

t  =  0  ; 

for  (i=0;  i<nsize  ;  i++) 

{ 

temp[i]  =  0; 

} 

while  (1)  //  Beginning  of  infinite  loop 

{ 

t++; 

//  generate  potential 

potent ial_generat ion (t)  ; 

for  (i=0;  i<nsize  ;  i++) 

{ 

vnoise[i]  =  v[i]  ; 

} 

//  add  noise 

add_noise (snr , idum)  ; 

//  For  1st  time  step  set  prev.  potential  to  cur.  pot. 

if  (t  :==  1) 

{ 

for  {i=0;  i<nsize;  i++) 

{ 

v0[i]  =  vnoise[i]  ; 

} 

) 

//  call  Hilde  filter 

flag  =  0  ; 
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if  (t  ==  1)  flag  =  1  ; 
hilde (vnoise, vO , hil, f lag)  ; 
rect (hil, hrec)  ; 
for  (i=0;  i<nsize;  i++) 

{ 

if  (potopt) 

{ 

potr[i]  =  hrec[i]  ; 

} 

else 

{ 

potr[i]  =  vnoise[i]  ; 

} 

} 

//  time  stepping  for  nt  steps  before  updating  the  pot. 

for  {it=l  ;  it  <  nt+1  ;  it++) 

{ 

time_stepping(psir,psii, nsize^mass)  ; 

//  calculate  the  probability 

for  (i=0;  i<nsize  ;  i++) 

.  { 

prob[i]  =  psir [i] *psir [i]  +  psii [ i] *psii [ i]  ; 

} 


//  double  buffer  control 

if  (page) 

page  =  0; 

else 

page  =  1 ; 

setactivepage (page) ; 
delay (ide lay) ; 
cleardevice ( ) ; 

//  plot  results 

Draw  (x_target ,  y_target ,  v,  scale_target )  ; 

Draw  (x_target , y_target ,  temp,  scale_target )  ; 

Draw (x_noise ,  y_noise ,  vnoise ,  scale_noise )  ; 

Draw  ( x_hi Ide ,  y_hi Ide ,  hi  1 ,  sea le_hi Ide )  ; 

Draw(x_pot,y_pot,hrec,  scale_pot)  ; 
Draw(x_prob,y_prob,  V,  scale_target )  ; 

Draw  (x_prob, y_prob,  temp,  scale_prob)  ; 
Draw(x__prob,y_prob,prob,  scale_prob)  ; 

Draw  ( x_pr ob ,  y_pr ob ,  po t  r ,  s  ca  1  e_^o t )  ; 

setvisualpage (page) ; 

) 

//  end  "it”  loop 
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if  ( i wopt ) 

{ 

//  calculate  overlap  with  ground  state 

sum  =  0 ; 
fr  =  0; 
f  i  =  0  ; 

for  (i=0;  i<nsize;  i++) 

{ 

sum  =  sum  +  psir [ i] *psir [i]  +  psii [ i ] *psii [ i ]  ; 
fr  =  fr  +  ground_state [i] *psir [i]  ; 
fi  =  fi  +  ground^state [i] *psii [i]  ; 

} 

if  (output) 

{ 

fout  «  ( f r*fr+f i*f i )  <<  "  "  «  sum  «  endl  ; 

} 


//  check  for  keyboard  hit 

if  (kbhit ()) break; 

}  //  end  of  infinite  loop 


/*  clean  up  */ 

getch ( ) ; 
closegraph ( ) ; 


C++  PROGRAM 


#include  <iostream.h> 
#include  <fstream.h> 

# include  <math.h> 


//  global  variables 

const  int  nsize  =  128  ; 
float  gain  =80.  ; 
float  dt  =  .05  ; 
int  nt  =  2  0  ; 
updates  of 
long  seed  =  3411  ; 


//  For  hilde  filter 
//  Size  of  time  step 
//  Number  of  time  step  between 
//  the  potential 

//  Seed  for  random  number  generator 


/*  Whenever  the  mass 
check  that  the  time 

float  beta  =  0.01  ; 
float  length  =  12.8  ; 


is  changed  (make  smaller)  remember 
step  dt  is  adequate  !  !  !  !  1  1  !  !  */ 

//  Mass  parameter 
//  Length  of  grid 


to 
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//Options 


int 

nopt  = 

1 

/ 

// 

1  = 

1  potential/  0=2  potential 

int 

iopt  = 

0 

7 

// 

1  = 

square  potential. 

// 

0  = 

l/cosh^2  potential 

int 

pot opt 

= 

0  ; 

// 

If  = 

=  1  use  hilde  filter  , 

// 

if  ; 

5  0  use  square  (cosh)  well 

int 

wopt  = 

0 

7 

// 

1  = 

constant  psi ( t=0 ) , 

// 

0  = 

ground  state  (for  cosh  pot) 

int 

output 

= 

1  ; 

// 

Output  overlap  data:  1  =  yes 

int 

kspace 

= 

1  ; 

// 

Use 

k-space?  1  =  yes  ,  0  =  no 

//  Data  about  potentials 

float  vel  =  0  ; 
float  xcenter  =  55  ; 
float  vdepth  =1.0  ; 
float  width  =  4  ; 
float  vel2  =  0  ; 

float  xcenter2  =70  ; 
float  vdepth2  =  1.0  ; 
float  width2  =  4  ; 

float  wa  =  4  ; 
float  Va  =  0  ; 

float  noise_factor  =  0  ; 
float  snr  =  5  ; 

//  Hilda  filter  function 

/* 

The  Hilde  filter  has  the  following  steps: 

1.  Take  the  difference  between  this  potential  and  the 
last  potential 

2.  Convolution  with  a  cone  filter  (here  taken  to  be  8 
nearest  points) 

3.  Add  a  little  of  this  info  to  prev*  info  (last  loop) 

*/ 


//  Velocity  (in  pixel  per  frame) 

//  Center  of  potential  initially 
//  Depth  of  potential 

//  Width  of  potential  in  pixels 

//  Velocity  of  2nd  potential  (in 
//  pixel  per  frame) 

//  Center  of  2nd  potential  initially 

//  Depth  of  2nd  potential 

//  Width  of  2nd  potential  in  pixels 

//  Width  of  absorbing  potential 

//  Depth  of  absorbing  potential 

//  0  =  no  noise  /  1  =  yes 

//  Signal  to  noise 


void  hilde (float  *xl,  float  *x0,  float  *out,  float  *al,  float  *a2,  float 
*a3,  int  nn) 

{ 

float  cone [9]  =  {1,2, 3, 4, 5, 4, 3, 2,1}  ; 

//normalize  cone  filter 
float  Slim  ; 

int  nf  =  5  ;  //nf  =  (size  of  array  "cone**  +  l)/2 
sum  =  0  ; 
int  i; 

for  {i=0  ;  i<2*nf-l  ;  i++) 

{ 

sum  +=  cone[i]  ; 
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} 

for  {i=0  ; 

{ 


cone [i] 

} 


i<2*nf-l 
i-  sum  ; 


i++) 


//  Take  the  difference 

float  *diff  ; 

diff  =  new  float [nn]  ; 

for  (i=0  ;  i<nn  ;  i++) 

{ 

diff[i]  =  xl[i]  -  xO[i]  ; 

} 


//  Convolution  with  a  cone  filter 

int  j; 

for  (i=0  ;  i<nn  ;  i++) 

{ 

sum  =  0  ; 

for  (j=:-{nf-l);  j  <  nf  ;  j++) 

{ 

sum  +=  cone [ j+nf “1] *dif f [ (nn+i+ j ) %nn]  ; 

} 

out  [  i  ]  =  sum  ; 


} 


//  Step  no,  3.  Hf actor  is  fraction  of  past  info  to  keep 

float  Hfactor  =  0.8  ; 
for  (i=0  ;  i<nn  ;  i++) 

{ 

al[i]  =  Hfactor^al [ i]  +  ( 1-Hfactor) *out [ i]  ; 
a2[i]  =  Hfactor*a2 [i]  +  ( 1-Hfactor ) *al [ i]  ; 
a3[i]  =  Hfactor*a3 [ i]  +  ( 1-Hfactor ) *a2 [ i ]  ; 
out [i]  =  a3 [i]  ; 

x0[i]  =  xl[i]  ; 

} 

delete[]  diff; 


/* 

Function  to  rectified  the  output  from  a  Hilde  filter, 
also  convolve  it  with  a  gaussian  3  times. 

*/ 


void  rect( float  *in,  float  *out,  int  nn,  float  gain) 

{ 

float  gauss[9]  =  {54,129,242,352,399,352,242,129,54}  ; 
//normalize  gaussian  filter 
float  sum  ; 

int  nf  =  5  ;  //  nf  =  (size  of  array  "gauss"  +  l)/2 

sum  =  0  ; 
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int  i ; 

for  (i=0  ;  i<2*nf-l  ;  i++) 

{ 

sum  +=  gauss [i]  ; 

} 

for  (i=0  ;  i<2*nf-l  ;  i++) 

{ 

gauss [i]  /=  sum  ; 

} 

float  *tem  ; 

tern  =  new  float [nn]  ; 

float  *tem2  ; 

tem2  =  new  float [nn]  ; 

//  Rectification  and  amplification 

for  (i=0  ;  i<nn  ;  i++) 

{ 

if  (in[i]  <  0) 

tem[i]  =  gain* in [i]  ; 

else 

tem[i]  =  0  ; 


//  Convolution  with  a  gaussian  filter  3 

int  j,k; 

for  (k=:0  ;  k<3  ;  k++) 

f 

for  (i=0  ;  i<nn  ;  i++) 

{ 

sum  =  0 ; 

for  (j=-(nf-l);  j  <  nf  ;  j++) 

{ 

sum  +=  gauss [j+nf”l] *tem[ (nn+i+j ) %nn] 

} 

tem2[i]  =  sum  ; 

} 

for  (i=0  ;  i<nn  ;  i++) 

{ 

tem[ i]  =  tem2 [i] ; 

} 

} 

for  (i=0  ;  i<nn  ;  i++) 

{ 

out[i]  =  tem[i]  ; 

} 

delete []  tern; 
delete []  tem2; 


times 
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//  Random  number  generator  (Adapted  from  Numerical  Recipes.*, 
by  Press  et .  al.) 

float  rani  (long  Scidum) 

{ 

const  long  Ml  =  259200  ; 
const  long  lAl  =  7141  ; 
const  long  ICl  =  54773  ; 
const  float  RMl  =  (1.0 /Ml)  ; 
const  long  M2  =  134456  ; 
const  long  IA2  =  8121  ; 
const  long  IC2  =  28411  ; 
const  float  RM2  =  (1.0/M2)  ; 
const  long  M3  =  243000  ; 
const  long  I  A3  =  4561  ; 
const  long  IC3  =  51349  ; 

static  long  ixl,  ix2,  ix3  ; 
static  float  r[98]  ; 
float  temp  ; 
static  int  iff=0  ; 
int  j  ; 

if  (idum  <  0  I  I  iff  ==  0) 

{ 

iff  =  1  ; 

ixl  =  (ICl  -  (idum))  %  Ml  ; 

ixl  =  (IAl*ixl  +  ICl)  %  Ml  ; 

ix2  =  ixl  %  M2  ; 

ixl  =  (lAl^ixl  +  ICl)  %  Ml  ; 

ix3  =  ixl  %  M3  ; 

for  (j=l  ;  j<=97;  j++) 

{ 

ixl  =  (IAl*ixl  +  ICl)  %  Ml  ; 
ix2  =  (IA2*ix2  +  IC2)  %  M2  ; 
r[j]  =  (ixl  +  ix2*RM2)*RMl  ; 

} 

idum  =  1  ; 

} 

ixl  =  (IAl*ixl  +  ICl)  %  Ml  ; 

ix2  :=  (IA2*ix2  +  IC2)  %  M2  ; 

ix3  :=  (IA3*ix3  +  IC3)  %  M3  ; 

j  =  1  +  ( {97*ix3)/M3)  ; 

temp  =  r[j]  ; 

r[j]  =  {ixl+ix2*RM2) '*^RM1  ; 
return  temp; 


//  fft  routine  (Adapted  from  Numerical  Recipes...  by  Press  et . 
al.  ) 

/*  A  section  was  added  at  the  end  to  divide  the  result  by  1/N 
if  the  inverse  FFT  is  taken  *  / 
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void  fourl (float  *data,  int  nn,  int  isign) 

( 

int  n,  inmax,  m,  j,  istep,  i  ; 

double  wtemp,  wr,  wpr,  wpi,  wi,  theta  ; 

float  tempr,  tempi,  tern  ; 


n  =  nn  «  1  ;  //n  =  nn*2 

j  =  1  ; 

for  (i=l;  i<n;  i+=2) 

{ 

if  (j  >  i) 

{ 

tern  =  data[i-l]  ; 
data[i-l]  =  data[j“l]  ; 
data[j“l]  =  tern  ; 
tern  =  data[i]  ; 
data[i]  =  data[j]  ; 
data[j]  =  tern  ; 

} 

m  =  n  »  1  ;  //m  =  n/2 
while  (m>=2  &&  j  >  m) 

{ 

j  -=  m  ;  //j  =  j  -  m 

m  »=  1  ;  //m  =  m/2 

} 

j  +=  m  ;  //j  =  j+m 

} 

mmax  =  2  ; 
while  (n  >  mmax) 

{ 

istep  =  2*mmax  ; 

theta  =  6.28318530717959/ (isign *mmax)  ; 

wtemp  =  sin ( 0 . 5* theta)  ; 

wpr  =  -2 . 0*wtemp*wtemp  ; 

wpi  =  sin (theta)  ; 

wr  =  1 . 0  ; 

wi  =  0 . 0  ; 

for  (m=l;  m  <  mmax;  m  +=  2) 

{ 

for  (i=m;  i<=n;  i  +=  istep) 

{ 

j  =  i  +  mmax  ; 

tempr  =  wr*data[j"l]  -  wi*data[j] 
tempi  =  wr*data[j]  +  wi*data[j-l] 
data[j-l]  =  data[i-l]  -  tempr  ; 
data[j]  =  data[i]  -  tempi  ; 
data[i-l]  +=  tempr  ; 
data[i]  +=  tempi  ; 


wtemp  =  wr  ; 

wr  =  wr*wpr  -  wi*wpi  +  wr 
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wi  =  wi'^'wpr  +  wtemp*wpi  +  wi  ; 

} 

mmax  =  istep  ; 

} 

//if  isign  =  -1  multiply  by  1/nn 
if  (isign  ==  -1) 

{ 

for  (i=0;  i  <  2*nn;  i++) 

{ 

data[i]  =  data[i]/nn  ; 

} 

} 

} 

//  Function  to  multiply  a  2  component  vector  with  a  2x2  matrix 
S 

void  Smultiply (float  Scpsilr,  float  Scpsili,  float  Scpsi2r,  float  &psi2i) 

{ 

double  sll  =  l./sqrt(2,)  ; 

double  sl2  =  l./sqrt(2.)  ; 

double  s21  =  l./sqrt(2.)  ; 

double  s22  =  -l./sqrt(2.)  ; 

double  temlr,  temli,  tem2r,  tem2i  ; 

temlr  =  sll*psilr  +  sl2*psi2r  ; 
temli  =  sll*psili  +  sl2*psi2i  ; 
tem2r  =  s21*psilr  +  s22*psi2r  ; 
tem2i  =  s21*psili  +  s22*psi2i  ; 

psilr  =  (float) (temlr)  ; 
psili  =  (float) (temli)  ; 
psi2r  =  (float) (tem2r)  ; 
psi2i  =  (float) (tem2i)  ; 

} 

//  Time  propagation  function 

/* 

parameters  to  be  passed  to  this  functions  are; 
number  of  grid  points  =  nn 

real  part  of  the  wave  function  psi  =  psir[nn] 

im  part  of  the  wave  function  psi  =  psii[nn] 

mass  =  mass 
time  step  =  dt 

real  part  of  potential  =  potr[nn] 

im  part  of  potential  =  poti  [nn] 

length  of  grid  =  length 

option  to  calculate  in  k  space  or  real  space  =  kspace 
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void  time_stepping  ( int  nn,  float  "^psir,  float  *psii,  float  mass,  float 
dt,  float  *potr,  float  *poti,  float  length,  int  kspace) 

{ 

int  i ,  j  ; 

float  alpha,  temr,  temi,  dx,  tem2r,  tem2i  ; 

float  psilr,  psili,  psi2r,  psi2i  ; 

float  *dat  ; 

dat  =  new  float [2*nn]  ; 

//  Set  up  the  kinetic  energy  vector 

float  *k2m  ; 

k2m  =  new  float  [nn]  ; 

double  pi  =  3.141592653589  ; 

f loat  dk  ; 

dk  =  (float) (2*pi/length)  ; 
for  (i=0;  i<nsi2e  ;  i++) 

{ 

k2in[i]  =  dk*(i)  ; 
if  { i  >  (nsize/2 ) ) 

{ 

k2m[i]  =  dk*(i-nsize)  ; 

} 

k2m[i]  =  k2m[ i] *k2m[i] /mass  ; 


dx  =  length/nn  ; 

//  Applied  exp[“iVt/2]  to  the  wave  function 

for  (i=0;  i<nn  ;  i++) 

{ 

alpha  =  potr[i]*dt/2  ; 

temr  =  cos (alpha) *psir [i]  +  sin (alpha) *psii [ i]  ; 
temi  =  cos (alpha) *psii [ i]  -  sin (alpha) *psir [ i ]  ; 

psir[i]  =  temr*exp(-poti [i] *dt/2)  ; 

psii[i]  =  temi*exp(-poti [i] *dt/2)  ; 

) 

if  (kspace)  //  Option  to  use  the  k-spcae  or  R- space  method 

{ 

//  First  fft  to  k  space 

for  (i=0;  i<nn  ;  i++) 

{ 

j  =  2*(i)  ; 

dat [ j ]  =  psir [i]  ; 

dat[j+l]  =  psii[i]  ; 

} 

fourl (dat , nn, 1)  ; 

for  (i=0;  i  <nn  ;  i++) 

{ 
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j  =  2*{i)  ; 
psir[i]  =  dat[j]  ; 
psii[i]  =  dat[j+l]  ; 


//  Applied  exp [-i (k**2 ) t /2m]  to  the  wave  function 

for  (i=0;  i<nn  ;  i++) 

{ 

alpha  =  k2m[i]*dt/2  ; 

temr  =  cos (alpha) *psir [i]  +  sin (alpha) *psii [ i ]  ; 
temi  =  cos (alpha) *psii [i]  -  sin (alpha) *psir [ i ]  ; 

psir[i]  =  temr  ; 
psii[i]  =  temi  ; 

} 


//  fft  back  to  real  space 

for  (i=0;  i<nn  ;  i++) 

{ 

j  =  2*{i)  ; 
dat[j]  =  psir[i]  ; 
dat[j+l]  =  psii[i]  ; 

} 

f our 1 ( dat , nn , - 1 )  ; 
for  (i=0;  i  <nn  ;  i++) 

{ 

j  =  2*(i)  ; 
psir[i]  =  dat[j]  ; 
psii[i]  =  dat[j+l]  ; 

} 

} 

else 

{ 

//  Everything  is  done  in  real  space 

//  Operate  exp [iAdt / (4*m*dx*dx) ] 

for  (i=0;  i<nn  ;  i  +=  2) 

{ 

psilr  =  psir[i]  ; 
psili  =  psii[i]  ; 
psi2r  =  psir[i+l]  ; 
psi2i  =  psii[i+l]  ; 

Smultiply (psilr, psili, psi2r, psi2i)  ; 

alpha  =  dt/4/mass/dx/dx* (-2)  ; 

tem2r  =  cos (alpha) *psi2r  -  sin (alpha) *psi2i  ; 

tem2i  =  cos (alpha) *psi2i  +  sin (alpha) *psi2r  ; 

psi2r  =  tem2r  ; 

psi2i  =  tem2i  ; 

Smultiply (psilr, psili, psi2r,psi2i)  ; 
psir[i]  =  psilr  ; 
psii  [ i]  =  psili  ; 
psir[i+l]  =  psi2r  ; 
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psii[i+l]  =  psi2i  ; 

} 

//  Operate  exp [ iBdt / (2*m*dx*dx) ] 

psilr  =  psirfO]  ; 
psili  =  psii [ 0 ]  ; 
psi2r  =  psir[nn“l]  ; 
psi2i  =  psii[nn“l]  ; 

Smultiply (psilr, psili, psi2r,psi2i)  ; 

alpha  =  dt/2/mass/dx/dx* (-2)  ; 

tem2r  =  cos  (alpha)  *psi2r  ~  sin  (alpha) '*'psi2i  ; 

tem2i  =  cos (alpha) *psi2i  +  sin (alpha) *psi2r  ; 

psi2r  =  tem2r  ; 

psi2i  =  tem2i  ; 

Smultiply (psilr, psili, psi2r, psi2i)  ; 

psir[0]  =  psilr  ; 

psii  [0]  =  psili  ; 

psir[nn‘“l]  =  psi2r  ; 

psii[nn-l]  =  psi2i  ; 

for  (i=l;  i<nn-l  ;  i  +=  2) 

{ 

psilr  =  psir[i]  ; 
psili  =  psii[i]  ; 
psi2r  =  psir[i+l]  ; 
psi2i  =  psii[i+l]  ; 

Smultiply (psilr, psili, psi2r,psi2i)  ; 

alpha  =  dt/2/mass/dx/dx* ( -2)  ; 

tem2r  =  cos (alpha) *psi2r  -  sin (alpha) *psi2i 

tem2i  =  cos (alpha) *psi2i  +  sin (alpha) *psi2r 

psi2r  =  tem2r  ; 

psi2i  =  tem2i  ; 

Smultiply (psilr, psili, psi2r,psi2i)  ; 

psir[i]  =  psilr  ; 

psii[i]  =  psili  ; 

psir[i+l]  =  psi2r  ; 

psii[i+l]  =  psi2i  ; 

} 

//  Operate  exp [iAdt / {4*m*dx*dx) ] 

for  (i=0;  i<nn  ;  i  +=  2) 

{ 

psilr  =  psir[i]  ; 
psili  =  psii[i]  ; 
psi2r  =  psir[i+l]  ; 
psi2i  =  psii[i+l]  ; 

Smultiply (psilr,psili,psi2r,psi2i)  ; 
alpha  =  dt/4/mass/dx/dx* (-2)  ; 
tem2r  =  cos (alpha) *psi2r  -  sin (alpha) *psi2i 
tem2i  =  cos (alpha) *psi2i  +  sin (alpha) *psi2r 
psi2r  =  tem2r  ; 
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psi2i  =  tem2i  ; 

Smultiply (psilr,psili,psi2r, psi2i)  ; 

psir[i]  =  psilr  ; 

psii [i]  =  psili  ; 

psir[i+l]  =  psi2r  ; 

psii [ i+1]  =  psi2i  ; 

} 

} 


//  Applied  exp[-iVt/2]  to  the  wave  function  again 

for  (i=0;  i<nn  ;  i++) 

{ 

alpha  =  potr[i]*dt/2  ; 

temr  =  cos (alpha) *psir [i]  +  sin (alpha) *psii [ i]  ; 
temi  =  cos (alpha) *psii [i]  -  sin(alpha) *psir [i]  ; 
psir[i]  =  teinr*exp(-poti  [i]  *dt/2)  ; 
psii[i]  =  temi*exp(-poti [i] *dt/2)  ; 

} 


delete []  dat  ; 
delete  []  k2in  ; 


//  Function  to  generate  the  potential 

void  potential_generat  ion  ( long  t,  float  int  nn) 

{ 


double  xupper,  xlower,  x,  xupper2,  xlower2,  y,  x2,  y2  ; 
int  i  ; 

if  (nopt)  //  Option  for  1  or  2  wells 

{ 

if  (iopt)  //  Option  for  square  or  cosh  well 

{ 


//  Set  up  square  potential 

xlower  =  vel*t  -  width/2  +  xcenter  ; 
xupper  =  vel*t  +  width/2  +  xcenter  ; 
xlower  =  f mod ( (xlower +nn)  ,  nn)  ; 
xupper  =  fmod ( (xupper +nn)  ,  nn)  ; 
for  (i=0  ;  i  <  nn  ;  i++) 

{ 

v[i]  =  0  ; 
if  (xlower  <  xupper) 

{ 

if  (i  >=  xlower  &&  i  <=  xupper)  v[i]  =  -vdepth; 

} 

else 

{ 

if  (i  <=  xupper  | |  i  >=  xlower)  v[i]  =  -vdepth; 

) 

} 

} 

else 
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//  Set  up  l/cosh**2  potential 

for  (i=0;  i<nn;  i++) 

{ 

X  =  i  -  xcenter  -  vel'^'t  ; 
y  =  (cosh (x/width) )  ; 
v[i]  =  -vdepth/y/y  ; 

} 

} 

} 

else 

{ 

if  (iopt)  //  Option  for  square  or  cosh  well 

{ 

//  Square  wall 

xlower  =  vel*t  -  width/2  +  xcenter  ; 
xupper  =  vel*t  +  width/2  +  xcenter  ; 
xlower  =  fmod ( (xlower+nn)  ,  nn)  ; 
xupper  =  fmod ( (xupper +nn)  ,  nn)  ; 
xlower2  =  vel2*t  -  width2/2  +  xcenter2  ; 
xupper2  =  vel2*t  +  width2/2  +  xcenter2  ; 
xlower2  =  fmod ( (xlower2+nn)  ,  nn)  ; 
xupper 2  =  fmod ( (xupper 2 +nn)  ,  nn)  ; 
for  (i=0  ;  i  <  nn  ;  i++) 

{ 

v[i]  =0  ; 

//  Check  to  see  if  point  falls  into  1st  well 

if  (xlower  <  xupper) 

{ 

if  (i  >=  xlower  &&  i  <=  xupper)  v[i]  =  -vdepth; 

} 

else 

{ 

if  (i  <=  xupper  II  i  >=  xlower)  v[i]  =  -vdepth; 

} 

//  Check  to  see  if  point  falls  into  2nd  well 

if  (xlower 2  <  xupper 2) 

{ 

if  (i  >=  xlower2  ScSc  i  <=  xupper2)  v[i]  =  -vdepth2 

} 

else 

{ 

if  (i  <=  xupper2  | |  i  >=  xlower2)  v[i]  =  -vdepth2 

} 

} 

} 

else 

{ 

//  Set  up  l/cosh**2  potential 

for  (i=0;  i<nn;  i++) 

{ 
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X  =  i  ““  xcenter  -  vel*t  ; 

y  =  (cosh (x/width) )  ; 

x2  =:  i  -  xcenter2  -  vel2*t  ; 

y2  =  (cosh (x2/width2 ) )  ; 

v[i]  =  -  vdepth/y/y  -  vdepth2/y2/y2  ; 

} 

} 

} 

} 

//  Main  program 

void  main  (void) 

{ 

int  i,  it  ; 
long  t ; 

float  sum,  mass,  factl,  fact2,  fr,  fi,  res,  x,  y  ; 

float  temp[nsize],  psir[nsize],  psii[nsize],  prob[nsi2e]  ; 

float  al[nsize]  ,  a2[nsize],  a3[nsize],  ground_state [nsize]  ; 

float  potr [nsize],  poti [nsize],  v[nsize] ,  vO [nsize],  vnoise [nsize] 

float  hil [nsize],  hrec [nsize]  ; 

double  z,  ss,  yy  ; 

long  idum  ; 

double  pi  =  3.141592653589  ; 
idum  =  seed  ; 

of St ream  f out ( "overlap . txt " )  ; 
of stream  vout { "potential . txt ” )  ; 

//  Calculate  the  mass  from  beta  (a  ratio  of  PE  to  KE) 

mass  =  beta* (pi^nsize/length) * (pi^nsize/length) /2/vdepth  ; 

if  (wopt) 

{ 

//  Initialize  wave  function  to  a  constant 

for  (i=0  ;  i<nsi2e  ;  i++) 

{ 

psir[i]  =  1/sqrt ( length)  ; 
psii[i]  =  0  ; 

} 

} 

else 

{ 

//  Ground  state  of  l/cosh**2  potential 

for  (i=0;  i<nsize  ;  i++) 

{ 

z  =  tanh (( i-xcenter) /width)  ; 
res  =  nsize/length  ; 

ss  =  0 . 5* ( -1+sqrt ( l+8*mass*vdepth*width*width/res/res ) )  ; 

yy  =  1  -  z*z  ; 

psir[i]  =  { float ) (pow(yy, ss/2 .) )  ; 
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psii[i]  =  0  ; 

} 

//  Normalize  wave  function 

sum  =  0 ; 

for  (i=0;  i<nsize  ;  i++) 

{ 

sum  =  sum  +  psir [i ] *psir [ i]  ; 

} 

sum  =  sqrt(sum)  ; 
for  (i=0;  i<nsize;  i++) 

{ 

psir[i]  =  psir [i] /sum  ; 
ground_state [i]  =  psir[i]  ; 

} 

} 


//  Set  up  the  absorbing  potential  at  the  end  of  the  grid 

for  {i=0;  i<nsize  ;  i++) 

{ 

factl  =  i*i/wa/wa  ; 

fact2  =  ( i~nsize+l) * (i-nsize+1) /wa/wa  ; 
if  (factl  >  30)  factl  =  30  ; 
if  (fact2  >  30)  fact2  =  30  ; 

poti[i]  =  Va  *  (exp (-factl)  +  exp(-fact2))  ; 

} 

t  =  0  ; 

for  (i=0;  i<nsize  ;  i++) 

{ 

temp[i]  =  0; 
al[i]  =  0  ; 
a2[i]  =0  ; 
a3[i]  =:  0  ; 
v[i]  -  0; 

} 

for  (t=l  ;  t<  100  ;  t++) 

{ 

/ / t++ ; 

//  Generate  potential 

potential_genGration (t , V, nsize)  ; 

for  (i=0;  i<nsize  ;  i++) 

{ 

vnoisG[i]  =  v[i]  ; 

} 


//  Add  noise 

X  =  (float) (pow(10 snr/10 .) )  ; 
y  =  vdepth/x  ; 
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for  (i=0  ;  i<nsize  ;  i++) 

{ 

vnoise[i]  =  v[i]  +  (rani ( idum) -0 . 5 ) *y*noise_f actor  ; 

} 

//  For  1st  time  step  set  previous  potential  to  present 
//  potential 

if  (t  ==  1) 

{ 

for  (i=0;  i<nsi2G;  i++) 

{ 

v0[i]  =  vnoise[i]  ; 

} 

} 

//  Call  Hilde  filter 

hilde ( vnoisG , vO , hil , al , a2 , a3 , nsize )  ; 

//  Rectify,  amplify,  and  convolve  output  of  Hilde  filter 

rect (hil, hrec,nsize, gain)  ; 

if  (potopt)  //  Option  to  choose  a  po.  for  the  Schrod.  eq. 

{ 

for  (i=0;  i<nsize;  i++) 

{ 

potr[i]  =  hrec[i]  ; 

} 

} 

else 

{ 

for  (i=0;  i<nsize;  i++) 

{ 

potr[i]  =  vnoise[i]  ; 

} 


//  Time  stepping  for  nt  steps  before  updating  the  pot. 

for  (it=l  ;  it  <  nt+1  ;  it++) 

{ 

tiine_stepping  (nsize,psir,psii,inass,  dt,potr,poti,  length,  kspace) 

//  Calculate  the  probability 

for  (i=0;  i<nsi2e  ;  i++) 

{ 

prob[i]  =  psir [i] *psir [i]  +  psii [i] *psii  [i]  ; 

} 

} 

//  End  **it"  loop 

if  (I wopt ) 

{ 

//  Calculate  overlap  with  ground  state 

sum  =  0; 
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fr  =  0; 
fi  =  0; 

for  (i=0;  i<nsize;  i++) 

{ 

sum  =  sum  +  psir [ i] *psir [ i]  +  psii [i] *psii [ i]  ; 
fr  =  fr  +  ground_state [i] *psir [i]  ; 
fi  =  fi  +  ground_state [i] *psii [i]  ; 

} 

if  (output) 

{ 

fout  «  (fr*fr+f i*f i)  «  “  “  «  sum  «  endl  ; 

} 

} 

}  //  End  of  t  loop 

//  Output  the  real  part  of  the  potential 

for  (i  =  0  ;  i<nsize  ;  i++) 

{ 

vout  «  i  «  "  "  «  potr[i]  «  endl  ; 

} 
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